Calibration Update

library(gotmtools)
## Loading required package: rLakeAnalyzer
## Warning: replacing previous import 'stats::filter' by 'dplyr::filter' when
## loading 'gotmtools'
library(LakeEnsemblR)
## 
## 
##   _          _        _____                          _     _ ____  
##  | |    __ _| | _____| ____|_ __  ___  ___ _ __ ___ | |__ | |  _ \ 
##  | |   / _` | |/ / _ |  _| | '_ \/ __|/ _ | '_ ` _ \| '_ \| | |_) |
##  | |__| (_| |   |  __| |___| | | \__ |  __| | | | | | |_) | |  _ < 
##  |_____\__,_|_|\_\___|_____|_| |_|___/\___|_| |_| |_|_.__/|_|_| \_\
##                                                                    
## 
##               https://github.com/aemon-j/LakeEnsemblR
## LakeEnsemblR version 1.0.9 (2021-08-26) is loaded
## WARNING: Your current version of LakeEnsemblR (v1.0.9) is ahead of the master branch version (1.0.5)
## Development version may have an unexpected behaviour
## 
## Attaching package: 'LakeEnsemblR'
## The following object is masked from 'package:gotmtools':
## 
##     analyse_strat
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(rLakeAnalyzer)
library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:dplyr':
## 
##     rename
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following objects are masked from 'package:reshape':
## 
##     colsplit, melt, recast
library(RColorBrewer)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:reshape':
## 
##     stamp
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(Metrics)
## 
## Attaching package: 'Metrics'
## The following object is masked from 'package:gotmtools':
## 
##     rmse
setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/")

ncdf <- "output/ensemble_output_all_models_31Aug21.nc"

model <- c("FLake", "GLM", "GOTM", "Simstrat", "MyLake")

spin_up <- 190



plot_heatmap(ncdf, model = model) +
  scale_colour_gradientn(limits = c(0, 32),
                         colours = rev(RColorBrewer::brewer.pal(11, "Spectral"))) + theme_classic()
## Extracted temp from output/ensemble_output_all_models_31Aug21.nc
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

plot_ensemble(ncdf, model = model, var = "ice_height")

fit <- calc_fit(ncdf, model = model, spin_up = spin_up)
fit
## $FLake
##       rmse       nse         r       bias      mae      nmae
## 1 2.324808 0.9060516 0.9571831 -0.7640633 1.477884 0.3071051
## 
## $GLM
##       rmse       nse         r      bias      mae      nmae
## 1 1.714442 0.9470689 0.9736236 0.1523328 1.136328 0.3467074
## 
## $GOTM
##      rmse       nse         r       bias     mae      nmae
## 1 1.85635 0.9379438 0.9700531 -0.3822369 1.26915 0.3418291
## 
## $Simstrat
##      rmse       nse         r      bias      mae      nmae
## 1 1.50802 0.9590475 0.9812579 0.1904871 1.161026 0.2203282
## 
## $MyLake
##       rmse       nse         r       bias       mae      nmae
## 1 1.409933 0.9642017 0.9820286 0.09948835 0.9535173 0.3288752
## 
## $ensemble_mean
##      rmse       nse         r       bias       mae      nmae
## 1 1.17144 0.9752881 0.9878914 -0.1224119 0.7751371 0.2228305
# out <- analyze_ncdf(ncdf, model, spin_up = 190)
# out$stats 

## Plot residuals
plist <- plot_resid(ncdf = "output/ensemble_output_all_models_14Aug21.nc", var = "temp")
## Extracted temp from output/ensemble_output_all_models_14Aug21.nc
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
ggarrange(plotlist = plist)

Schmidt Stability

  • Not much to say about Schmidt Stability, observations seem to track well with the mean at the very least
  • RMSE shows that all models track relatively similarly in terms of GOF with the exception of FLake which has an excessively high RMSE at 239.7
  • Flake Schmidt stability looks quite low after changing bathymetry to include the surface value only. If possible a double check on my code for this section would be much appreciated to ensure a mistake wasn’t made in the process.
ncdf <- "~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_31Aug21.nc"


######################## Calculating Schmidt Stability ################################


out <- load_var(ncdf = ncdf, var = "temp")
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
# out <- as.data.frame(out[[1]])
bathy <- read.csv('~/Dropbox/sunapee_LER_projections/LER_inputs/sunapee_hypso.csv')
colnames(bathy) <- c("depths", "areas")
ts.sch <- lapply(out, function(x) {
  ts.schmidt.stability(x, bathy = bathy, na.rm = TRUE)
})


## Reshape to data.frame
df <- melt(ts.sch, id.vars = 1)
colnames(df)[4] <- "model"
df$yday <- yday(df$datetime)
df$year <- year(df$datetime)
df <- filter(df, model != "FLake")


## CALCULATING FLAKE SCHMIDT STABILITY
listflake <- list(out[["FLake"]])
flake_depth <- glmtools::get_nml_value(nml_file = "~/Dropbox/sunapee_LER_projections/LER_calibration/FLake/flake.nml", arg_name = "depth_w_lk")
bathy <- filter(bathy, depths <= flake_depth)
ts.sch <- lapply(listflake, function(x) {
  ts.schmidt.stability(x, bathy = bathy, na.rm = TRUE)
})

dfflake <- melt(ts.sch, id.vars = 1)
colnames(dfflake)[4] <- "model"
dfflake$model <- "FLake"
dfflake$yday <- yday(dfflake$datetime)
dfflake$year <- year(dfflake$datetime)

colnames(df)
## [1] "datetime" "variable" "value"    "model"    "yday"     "year"
colnames(dfflake)
## [1] "datetime" "variable" "value"    "model"    "yday"     "year"
df <- rbind(df, dfflake)
########################## END CALCULATING FLAKE SCHMIDT STABILITY


wideform <- dcast(df, datetime~model, value.var = "value")
wideform <- filter(wideform, is.na(Obs) == FALSE & is.na(GLM) == FALSE &
                     is.na(GOTM) == FALSE & is.na(FLake) == FALSE & 
                     is.na(Simstrat) == FALSE & is.na(MyLake) == FALSE)


rmse <- c()
models <- c("FLake", "GOTM", "Simstrat", "MyLake", "GLM")
rmse <- c(rmse, (rmse(wideform$Obs, wideform$FLake)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$GOTM)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$Simstrat)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$MyLake)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$GLM)))

models <- as.data.frame(models)
rmse <- as.data.frame(rmse)

model_rmse <- cbind(models, rmse)
model_rmse
##     models      rmse
## 1    FLake 166.16345
## 2     GOTM  69.78399
## 3 Simstrat  75.21572
## 4   MyLake  55.53663
## 5      GLM 113.50085
dfobs <- filter(df, model == "Obs")
df <- filter(df, model != "Obs")


df <- df %>% 
  dplyr :: group_by(yday, year) %>% 
  dplyr :: mutate(mean = mean(value, na.rm = TRUE)) %>% 
  dplyr :: mutate(sd = sd(value, na.rm = TRUE))

dfobs$mean <- NA
dfobs$sd <- NA

df <- rbind(df, dfobs)

## plot results
ggplot(df, aes(yday, value, colour = model)) +
  # geom_line(data=df, aes(y=mean, x=yday), color = "black", size = 1) +
  geom_line() +
  # geom_ribbon(data = df, aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6, 
  #             linetype = 0.1, 
  #             color = "grey") +
  facet_wrap(~year) +
  labs(y = "Schmidt stability (J/m2)") +
  theme_classic() + ylim(-50, 1000)

ggplot(subset(df, model != "Obs"), aes(yday, mean)) + 
  geom_line() + 
  geom_ribbon(data = subset(df, model != "Obs"), aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6,
              linetype = 0.1,
              color = "grey") + 
  facet_wrap(~year) + 
  labs(y = "Schmidt stability (J/m2)") +
  theme_classic() + ylim(-50, 1000) + 
  geom_line(data = subset(df, model == "Obs"), aes(yday, value, col = "Obs"))

Thermocline Depth

  • Similar update on thermocline depth.
  • GLM appears to be the poorest performer most years based on the first figure.
  • The extra years of observations give me more confidence, as 2009 and 2010 were a bit off. The thermocline seems to fit quite well from 2011 and on.
  • FLake is the highest performer for RMSE, with MyLake being the lowest performer. When looking at MyLake on the plots however, it seems to perform relatively well with occassional anomalous jumps down to the max depth towards the end of the season. GLM may perform worse than MyLake during the overall season as the RMSE is 7 and GLM consistently models significantly lower than the observed thermocline depth.
library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

ncdf <- "~/Dropbox/sunapee_LER_projections/LER_calibration/output/ensemble_output_all_models_31Aug21.nc"

out <- load_var(ncdf = ncdf, var = "temp")
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
## Same for thermocline depth
ts.td <- lapply(out, function(x) {
  ts.thermo.depth(x, Smin = 0.1, na.rm = TRUE)
})

df <- melt(ts.td, id.vars = 1)
colnames(df)[4] <- "model"
df$yday <- yday(df$datetime)
df$month <- month(df$datetime)
df$year <- year(df$datetime)

wideform <- dcast(df, datetime~model, value.var = "value")
wideform <- filter(wideform, is.na(Obs) == FALSE & is.na(GLM) == FALSE &
                     is.na(GOTM) == FALSE & is.na(FLake) == FALSE & 
                     is.na(Simstrat) == FALSE & is.na(MyLake) == FALSE)


rmse <- c()
models <- c("FLake", "GOTM", "Simstrat", "MyLake", "GLM")
rmse <- c(rmse, (rmse(wideform$Obs, wideform$FLake)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$GOTM)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$Simstrat)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$MyLake)))
rmse <- c(rmse, (rmse(wideform$Obs, wideform$GLM)))

models <- as.data.frame(models)
rmse <- as.data.frame(rmse)

model_rmse <- cbind(models, rmse)
model_rmse
##     models      rmse
## 1    FLake  2.283110
## 2     GOTM  4.573914
## 3 Simstrat  3.415918
## 4   MyLake 13.902624
## 5      GLM  6.597361
dfobs <- filter(df, model == "Obs")
df <- filter(df, model != "Obs")

df <- df %>% 
  dplyr :: group_by(yday, year) %>% 
  dplyr :: mutate(mean = mean(value, na.rm = TRUE)) %>% 
  dplyr :: mutate(sd = sd(value, na.rm = TRUE))

dfobs$mean <- NA
dfobs$sd <- NA

df <- rbind(df, dfobs)


ggplot(subset(df, month >= 6 & month <= 8), aes(yday, value, colour = model)) +
  facet_wrap(~year) +
  geom_line() +
  labs(y = "Thermocline depth (m)") +
  scale_y_continuous(trans = "reverse") +
  theme_classic() 

ggplot(subset(df, month >= 6 & month <= 8 & model != "Obs"), aes(yday, mean)) +
  facet_wrap(~year) +
  geom_line() +
  geom_ribbon(data = subset(df, month >= 6 & month <= 8 & model != "Obs"), aes(ymin = mean-sd, ymax = mean+sd), alpha = 0.6,
              linetype = 0.1,
              color = "grey") + 
  labs(y = "Thermocline depth (m)") +
  scale_y_continuous(trans = "reverse") +
  theme_classic() + 
  geom_line(data = subset(df, month >= 6 & month <= 8 & model == "Obs"), aes(yday, value, col = "Obs"))

dfsummer <- filter(df, month >=6 & month <=8)

All other metrics

  • Observations are added for this group of metrics
  • There seems to be an issue with 2009 calculations, as shown by the printed dataframe.
  • The only metric to be compared with ice will have to be added in - coming soon.
  • There is very limited stratification data calculated for the observed data.
  • Seems to very high uncertainty for bottom temperature. Reason to remove it from desired metrics?
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
##            [,1]                  
## short_name "ice_height"          
## units      "m"                   
## dimensions "lon lat model member"

## # A tibble: 24 x 6
## # Groups:   year, variable [24]
##     year variable      value model  mean    sd
##    <dbl> <fct>         <dbl> <chr> <dbl> <dbl>
##  1  2009 TsMax         0.140 Obs      NA    NA
##  2  2009 TsMaxDay      5     Obs      NA    NA
##  3  2009 TsMin        -0.270 Obs      NA    NA
##  4  2009 TsMinDay      0     Obs      NA    NA
##  5  2009 TbMax         0.970 Obs      NA    NA
##  6  2009 TbMaxDay     11     Obs      NA    NA
##  7  2009 TbMin         0.400 Obs      NA    NA
##  8  2009 TbMinDay      0     Obs      NA    NA
##  9  2009 MaxStratDur  NA     Obs      NA    NA
## 10  2009 MeanStratDur NA     Obs      NA    NA
## # … with 14 more rows

Ice off

  • GLM does not model ice so is not present
  • FLake, MyLake and Simstrat do extremely well when comparing by plot and by RMSE, with all models hovering slightly above or below 5 RMSE.
  • GOTM is performing poorly when modeling ice, with an RMSE of 59 and a consistently lower ice-off date than any of the other models.
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
##            [,1]                    
## short_name "temp"                  
## units      "degC"                  
## dimensions "lon lat z model member"
##            [,1]                  
## short_name "ice_height"          
## units      "m"                   
## dimensions "lon lat model member"
##     models      rmse
## 1    FLake  4.459696
## 2     GOTM 59.897134
## 3 Simstrat  4.333333
## 4   MyLake  4.887626

library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/calibration_results_FLake_083021_v2/")

out_f <- "calibration_results_FLake_083021_v2"
config_file <- 'LakeEnsemblRsun.yaml'
model <- c("FLake")

param_file <- "calibration_results_FLake_083021_v2/FLake_LHC_202108311031.csv"


cal_files <- list.files(".", full.names = TRUE)
cal_files <- cal_files[c(1,2)]

res <- load_LHC_results(config_file = config_file, model = model, res_files = cal_files)
dim(res[[model]])
## [1] 500  12
df <- plyr::ldply(res, function(x) {
  df <- x[, -c(3:7)]
  reshape2::melt(df, id.vars = c("par_id", "rmse"))
}, .id = "model")
df$id_no <- as.numeric(gsub(".*?([0-9]+).*", "\\1", df$par_id)) 

bst_par <- df$id_no[which.min(df$rmse)]
sub <- df[df$id_no == bst_par, ]
sub
##      model par_id     rmse                variable      value id_no
## 330  FLake  p0330 2.324808              wind_speed  1.0404000   330
## 830  FLake  p0330 2.324808                     swr  0.8437100   330
## 1330 FLake  p0330 2.324808   LAKE_PARAMS/c_relax_C  0.0054947   330
## 1830 FLake  p0330 2.324808 LAKE_PARAMS/depth_bs_lk  2.4575000   330
## 2330 FLake  p0330 2.324808     LAKE_PARAMS/T_bs_lk 15.2050000   330
p1 <- ggplot(df) +
  geom_point(aes(value, rmse)) +
  facet_wrap(model~variable, scales = "free_x") +
  geom_hline(yintercept = 2, linetype = "dashed") +
  ylab("RMSE (\u00B0C)") +
  geom_vline(data = sub, aes(xintercept = value)) +
  geom_hline(yintercept = 3.5, color = "red", linetype = "dashed") +
  # coord_cartesian(ylim = c(1, 4)) +
  # scale_x_log10() +
  theme_classic(base_size = 16)
p1

library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/calibration_results_GLM_082621/")

out_f <- "calibration_results_GLM_082621"
config_file <- 'LakeEnsemblRsun.yaml'
model <- c("GLM")

# param_file <- "calibration_results_FLake_083021_v2/FLake_LHC_202108311031.csv"


cal_files <- list.files(".", full.names = TRUE)
cal_files <- cal_files[c(1,2)]

res <- load_LHC_results(config_file = config_file, model = model, res_files = cal_files)
dim(res[[model]])
## [1] 500  12
df <- plyr::ldply(res, function(x) {
  df <- x[, -c(3:7)]
  reshape2::melt(df, id.vars = c("par_id", "rmse"))
}, .id = "model")
df$id_no <- as.numeric(gsub(".*?([0-9]+).*", "\\1", df$par_id)) 

bst_par <- df$id_no[which.min(df$rmse)]
sub <- df[df$id_no == bst_par, ]
sub
##      model par_id     rmse                   variable    value id_no
## 36     GLM  p0036 1.714442                 wind_speed  1.49590    36
## 536    GLM  p0036 1.714442                        swr  1.06550    36
## 1036   GLM  p0036 1.714442 sediment/sed_temp_mean...3  3.48260    36
## 1536   GLM  p0036 1.714442 sediment/sed_temp_mean...4 10.59700    36
## 2036   GLM  p0036 1.714442  glm_setup/max_layer_thick  0.64125    36
p1 <- ggplot(df) +
  geom_point(aes(value, rmse)) +
  facet_wrap(model~variable, scales = "free_x") +
  geom_hline(yintercept = 2, linetype = "dashed") +
  ylab("RMSE (\u00B0C)") +
  geom_vline(data = sub, aes(xintercept = value)) +
  geom_hline(yintercept = 3.5, color = "red", linetype = "dashed") +
  # coord_cartesian(ylim = c(1, 4)) +
  # scale_x_log10() +
  theme_classic(base_size = 16)
p1
## Warning: Removed 5 rows containing missing values (geom_point).

library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/calibration_results_GOTM_081321/")

out_f <- "calibration_results_GOTM_081321"
config_file <- 'LakeEnsemblRsun.yaml'
model <- c("GOTM")

# param_file <- "calibration_results_FLake_083021_v2/FLake_LHC_202108311031.csv"


cal_files <- list.files(".", full.names = TRUE)
cal_files <- cal_files[c(1,2)]

res <- load_LHC_results(config_file = config_file, model = model, res_files = cal_files)
dim(res[[model]])
## [1] 500  10
df <- plyr::ldply(res, function(x) {
  df <- x[, -c(3:7)]
  reshape2::melt(df, id.vars = c("par_id", "rmse"))
}, .id = "model")
df$id_no <- as.numeric(gsub(".*?([0-9]+).*", "\\1", df$par_id)) 

bst_par <- df$id_no[which.min(df$rmse)]
sub <- df[df$id_no == bst_par, ]
sub
##      model par_id    rmse                    variable      value id_no
## 56    GOTM  p0056 1.85635                  wind_speed 5.4340e-01    56
## 556   GOTM  p0056 1.85635                         swr 1.3701e+00    56
## 1056  GOTM  p0056 1.85635 turbulence/turb_param/k_min 2.6656e-06    56
p1 <- ggplot(df) +
  geom_point(aes(value, rmse)) +
  facet_wrap(model~variable, scales = "free_x") +
  geom_hline(yintercept = 2, linetype = "dashed") +
  ylab("RMSE (\u00B0C)") +
  geom_vline(data = sub, aes(xintercept = value)) +
  geom_hline(yintercept = 3.5, color = "red", linetype = "dashed") +
  # coord_cartesian(ylim = c(1, 4)) +
  # scale_x_log10() +
  theme_classic(base_size = 16)
p1

library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/calibration_results_MyLake_081321/")

out_f <- "calibration_results_MyLake_081321"
config_file <- 'LakeEnsemblRsun.yaml'
model <- c("MyLake")

# param_file <- "calibration_results_FLake_083021_v2/FLake_LHC_202108311031.csv"


cal_files <- list.files(".", full.names = TRUE)
cal_files <- cal_files[c(1,2)]

res <- load_LHC_results(config_file = config_file, model = model, res_files = cal_files)
dim(res[[model]])
## [1] 500  10
df <- plyr::ldply(res, function(x) {
  df <- x[, -c(3:7)]
  reshape2::melt(df, id.vars = c("par_id", "rmse"))
}, .id = "model")
df$id_no <- as.numeric(gsub(".*?([0-9]+).*", "\\1", df$par_id)) 

bst_par <- df$id_no[which.min(df$rmse)]
sub <- df[df$id_no == bst_par, ]
sub
##       model par_id     rmse           variable   value id_no
## 219  MyLake  p0219 1.409933         wind_speed 1.48920   219
## 719  MyLake  p0219 1.409933                swr 1.27410   219
## 1219 MyLake  p0219 1.409933 Phys.par/C_shelter 0.39897   219
p1 <- ggplot(df) +
  geom_point(aes(value, rmse)) +
  facet_wrap(model~variable, scales = "free_x") +
  geom_hline(yintercept = 2, linetype = "dashed") +
  ylab("RMSE (\u00B0C)") +
  geom_vline(data = sub, aes(xintercept = value)) +
  geom_hline(yintercept = 3.5, color = "red", linetype = "dashed") +
  # coord_cartesian(ylim = c(1, 4)) +
  # scale_x_log10() +
  theme_classic(base_size = 16)
p1

library(gotmtools)
library(LakeEnsemblR)
library(ggplot2)
library(LakeEnsemblR)
library(ggpubr)
library(dplyr)
library(rLakeAnalyzer)
library(reshape)
library(RColorBrewer)
library(lubridate)

setwd("~/Dropbox/sunapee_LER_projections/LER_calibration/calibration_results_Simstrat_081321/")

out_f <- "calibration_results_Simstrat_081321"
config_file <- 'LakeEnsemblRsun.yaml'
model <- c("Simstrat")

# param_file <- "calibration_results_FLake_083021_v2/FLake_LHC_202108311031.csv"


cal_files <- list.files(".", full.names = TRUE)
cal_files <- cal_files[c(1,2)]

res <- load_LHC_results(config_file = config_file, model = model, res_files = cal_files)
dim(res[[model]])
## [1] 500  10
df <- plyr::ldply(res, function(x) {
  df <- x[, -c(3:7)]
  reshape2::melt(df, id.vars = c("par_id", "rmse"))
}, .id = "model")
df$id_no <- as.numeric(gsub(".*?([0-9]+).*", "\\1", df$par_id)) 

bst_par <- df$id_no[which.min(df$rmse)]
sub <- df[df$id_no == bst_par, ]
sub
##         model par_id     rmse                 variable     value id_no
## 20   Simstrat  p0020 1.502343               wind_speed 1.4883000    20
## 520  Simstrat  p0020 1.502343                      swr 1.0299000    20
## 1020 Simstrat  p0020 1.502343 ModelParameters/a_seiche 0.0035929    20
p1 <- ggplot(df) +
  geom_point(aes(value, rmse)) +
  facet_wrap(model~variable, scales = "free_x") +
  geom_hline(yintercept = 2, linetype = "dashed") +
  ylab("RMSE (\u00B0C)") +
  geom_vline(data = sub, aes(xintercept = value)) +
  geom_hline(yintercept = 3.5, color = "red", linetype = "dashed") +
  # coord_cartesian(ylim = c(1, 4)) +
  # scale_x_log10() +
  theme_classic(base_size = 16)
p1